source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

inflam15 <- read_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

# CRC -------------
sobj <- read_rds('CRC-I/data/crc_merge4_immune.rds')

sobj |> DimPlot(cols = DiscretePalette(36))

sobj |> DotPlot(pbmc_markers, cols = 'RdBu', cluster.idents = T)

sobj |> DimPlot(group.by = 'manual.main')

sobj |> DotPlot2d('TRPM2', manual.main, genotype) +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'FCGR2B-I232T')

sobj |> 
  mutate(dataset = str_to_title(dataset),
         manual.main = case_when(manual.main != 'Myeloid_cells' ~ manual.main,
                                 .cell %in% colnames(sobj_mf) ~ 'Macrophages',
                                 .default = 'DC')) |>
  DotPlot2d('TRPM2', manual.main, dataset) +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Dataset')

## filter myeloid ----------
sobj_myl <- sobj |>
  filter(manual.main == 'Myeloid_cells')

sobj_myl <- sobj_myl |>
  AddModuleScore(features = list(inflam15$gene), name = 'inflam')

sobj_myl <- sobj_myl |>
  quick_process_seurat(batch = c('orig.ident', 'dataset'), skip_norm = T)

modc_marker <- list('cMono'=c('CD14','S100A8','S100A9'),
                    'ncMono' = c('FCGR3A','CDKN1C'),
                    'Macrophage'=c('CD68','MAFB','MARCO','CD163','MRC1'),
                    'cDC'=c('FCER1A','CLEC9A','XCR1','CD1C','CLEC10A'),
                    'pDC'='TLR7')

sobj_myl |>
  DotPlot(c(modc_marker, 'TRPM2'), cols = 'RdBu') +
  RotatedAxis()

sobj_myl <- sobj_myl |>
  mutate(manual_fine = case_when(seurat_clusters %in% c(14,6,12) ~ 'cDC',
                                 seurat_clusters %in% c(13,11) ~ 'pDC',
                                 seurat_clusters %in% c(9) ~ 'CD4',
                                 .default = 'Macro'))

sobj_myl |>
  filter(manual_fine != 'CD4') |>
  mutate(dataset = str_to_title(dataset)) |>
  DotPlot2d('TRPM2', manual_fine, dataset) +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Dataset')

sobj_myl |>
  filter(manual_fine != 'CD4') |>
  DimPlot(group.by = 'manual_fine') +
  ggtitle('CRC TME myeloid cells')

sobj_myl |>
  filter(manual_fine != 'CD4') |>
  FeaturePlot('TRPM2', order = T, cols = c('lightgrey','red'))

sobj_mf <- sobj_myl |>
  filter(manual_fine == 'Macro')

sobj_mf <- sobj_mf |>
  FindClusters(algorithm = 4, random.seed = 1)

### by dataset ----------
m2_inflam <- sobj_mf |>
  filter(dataset == 'li2021') |>
  DotPlot(c('TRPM2','inflam1')) |>
  pluck('data') |>
  as_tibble()

m2_inflam |>
  filter(features.plot == 'TRPM2') |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  labs(title = 'CRC macrophage', subtitle = 'Li2021')

m2_inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .12, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2, label.x = .15) +
  labs(title = 'TME myeloid clusters in CRC tumor', subtitle = 'Li2021',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('CRC.mf.trpm2.inflam.correlation.li21', width = 70)

## inflam violin ---------
sobj_mf <- sobj_mf |>
  mutate(m2_type = ifelse(seurat_clusters == 12, 'TRPM2-hi MF', 'TRPM2-lo MF'))

sobj_mf |>
  DotPlot('TRPM2', cols = 'RdYlBu', dot.scale = 2) +
  labs(x = 'Gene', y = 'Macrophage clusters') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('CRC.mf.TRPM2.dotplot')

sobj_mf |>
  bill.violin(inflam15$gene, group.by = m2_type, facet.ncol = 5) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in CRC TME MF') +
  theme_pubr(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_pdf('CRC.m2hl.inflam.gene.violin.pdf', width = 90)

## immunosuppresive -----------
tam_immunosuppressive_signature <- c(
  # 精氨酸代谢
  "ARG1", "NOS2",
  # 色氨酸代谢
  "IDO1", "IDO2", "TDO2",
  # 免疫检查点配体
  "CD274", "PDCD1LG2", "LGALS9", "VSIR", "CD80", "CD86",
  # 抗炎/修复因子
  "IL10", "TGFB1", "IL4R", "IL1RN",
  # 趋化因子
  "CCL17", "CCL18", "CCL22", "CXCL12",
  # 吞噬抑制
  "SIRPA", "MERTK"
)

sobj_mf <- sobj_mf |>
  AddModuleScore(features = list(tam_immunosuppressive_signature),
                 name = 'immunosup')

m2_immunosup <- sobj_mf |>
  DotPlot(c('TRPM2','immunosup1')) |>
  pluck('data') |>
  as_tibble()

m2_immunosup |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .75, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  ggplot(aes(TRPM2, immunosup1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  #stat_cor(size = 2, label.x = .4, label.y = 1.4) +
  labs(title = 'TME myeloid clusters in CRC tumor',
       y = 'Immunosuppressive module score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('CRC.mf.trpm2.inflam.correlation', width = 70)

sobj_mf |>
  mutate(subtype = ifelse(seurat_clusters == 9, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  VlnPlot(c('immunosup1'), group.by = 'subtype', pt.size = 0) +
  labs(x = 'Subtype', title = 'Immunosuppresive score')

# ESCC -----------
sobj <- read_rds('CRC-I/esophagealCancer/zhang2021escc.rds')

sobj <- sobj |>
  mutate(type.main = case_when(str_detect(type.fine, 'B$|GCB|Plasma') ~ 'B',
                               str_detect(type.fine, 'DC') ~ 'DC',
                               str_detect(type.fine, 'Mono') ~ 'Mono',
                               str_detect(type.fine, 'TAM') ~ 'TAM',
                               type.fine == 'Mast' ~ 'Mast',
                               .default = 'T'))

sobj |> DimPlot(group.by = 'type.main') +
  ggtitle('ESCC TME')

sobj |> DotPlot2d('TRPM2', genotype, type.main) +
  labs(x = 'FCGR2B-I232T', y = 'Cell type')

m2_bysample <- sobj |> DotPlot2d('TRPM2', sample, type.main) |>
  pluck('data') |>
  as_tibble() 

m2_bysample |>
  mutate(group.y = fct_reorder(group.y, avg.exp)) |>
  ggplot(aes(group.y, avg.exp)) +
  geom_boxplot()

m2_bysample |>
  mutate(group.y = fct_reorder(group.y, avg.exp, .fun = mean)) |>
  ggplot(aes(group.y, avg.exp, fill = group.y)) +
  stat_summary(geom = 'col', fun = 'mean') +
  geom_jitter(height = 0, width = .1) +
  theme_bw() +
  NoLegend() +
  labs(x = 'Cell type', y = 'Average expression',
       title = 'TRPM2 expression in ESCC TME')

sobj_mf <- sobj |> filter(type.main == 'TAM')

sobj_mf <- sobj_mf |>
  quick_process_seurat(skip_norm = T)

sobj_mf <- sobj_mf |>
  AddModuleScore(features = list(inflam15$gene), name = 'inflam')

m2_inflam <- sobj_mf |>
  DotPlot(c('TRPM2','inflam1')) |>
  pluck('data') |>
  as_tibble()

m2_inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .2, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'TME macrophage clusters in ESCC',
       y = 'Inflammatory module score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('ESCC.mf.trpm2.inflam.correlation', width = 70)

# pan-cancer myeloid ------------
pcancer_supp <-
  GEOquery::getGEOSuppFiles('GSE154763', makeDirectory = F, fetch_files = F)

pcancer_supp |>
  as_tibble() |>
  filter(!str_detect(fname, 'cDC')) |>
  select(url) |>
  write_csv('mission/SLE_TRPM2_MfMo/data/z2021_pan_myl.lst', col_names = F)

library(data.table)

pmeta <- list.files('mission/SLE_TRPM2_MfMo/data/z2021_pan_myl', 'metadata',
                    full.names = T)

pexpr <- list.files('mission/SLE_TRPM2_MfMo/data/z2021_pan_myl', 'expr',
                    full.names = T)

logtpm2count <- function(norm_expr, meta, name){
  esca_meta <- meta |> read_csv(show_col_types = F) |>
    column_to_rownames('index')
  
  scale_factor <- median(esca_meta$n_counts)
  
  message(str_glue('Scale factor: {scale_factor}'))
  
  message('Reading log-TPM matrix...')
  esca_dense <- norm_expr |>
    fread() |>
    transpose(keep.names = 'gene', make.names = 'index') |>
    column_to_rownames('gene') |>
    as.sparse()
  
  esca_7676 <- expm1(esca_dense) / scale_factor
  
  cell_count_sum <- esca_meta$n_counts
  
  message('Value before rounding:')
  print(head(esca_7676@x))
  
  esca_round <- esca_7676 %*% Matrix::Diagonal(x = cell_count_sum)
  
  message('Recovered value:')
  print(head(esca_round@x))
  
  esca_round@x <- esca_round@x |> round()
  
  message('Final count value:')
  print(head(esca_round@x))

  colnames(esca_round) <- colnames(esca_7676)
  
  h5_path <- str_glue('{name}.h5')
  
  if(!file.exists(h5_path)) {
    esca_round |> DropletUtils::write10xCounts(path = str_glue('{name}.h5'))
  }
}

pp_add_module <- function(h5_path, meta, module_genes, module_name) {
  esca_meta <- meta |> read_csv(show_col_types = F) |>
    column_to_rownames('index')
  
  sobj <- h5_path |>
    Read10X_h5() |>
    CreateSeuratObject(meta.data = esca_meta) |>
    quick_process_seurat(batch = 'patient') |>
    AddModuleScore(features = list(module_genes), name = module_name)
  }

pname <- pmeta |>
  str_extract('(?<=763_).+(?=_meta)')

score1 <- logtpm2count(norm_expr = pexpr[1], meta = pmeta[1], 
                          name = pname[1])

score1

pan_cancer_m2inflam <- list(norm_expr = pexpr, meta = pmeta, name = pname) |>
  pmap(logtpm2count, .progress = T)

## macrophage ------------
ph5 <- list.files('mission/SLE_TRPM2_MfMo/data/z2021_pan_myl', 'h5',
                    full.names = T)

sobj_lst <- list(x = ph5, y = pmeta) |>
  pmap(\(x,y)pp_add_module(h5_path = x, meta = y,
                           module_genes = inflam15$gene, module_name = 'inflam'),
       .progress = T)

sobj_lst |>
  write_rds('mission/SLE_TRPM2_MfMo/data/z2021_pan_myl/pan_cancer_lst.rds')

sobj <- sobj_lst[[1]]

sobj_lst <- sobj_lst |>
  set_names(pname)

corr_in_celltype <- function(sobj, cell_regex, corr_gene = c('TRPM2','inflam1'),
                             redo_pca = TRUE){
  sobj <- sobj |>
    filter(str_detect(MajorCluster, cell_regex)) 
  
  if (redo_pca) {
    sobj <- sobj |>
      quick_process_seurat(skip_norm = T, batch = 'patient')
  } else {
    sobj <- sobj |>
      FindClusters(algorithm = 4, random.seed = 1)
  }
  sobj |>
    DotPlot(corr_gene) |>
    pluck('data') |>
    as_tibble()
  }

### TRPM2 x inflammatory ------------
pan_cancer_m2inf_mf <- sobj_lst |>
  map(\(x)corr_in_celltype(x, 'Macro', redo_pca = F), .progress = T)

pan_cancer_m2inf_mf |>
  list_rbind(names_to = 'cancer') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = c(cancer,id)) |>
  summarise(pearson = cor.test(TRPM2, inflam1) |> broom::tidy(), .by = cancer) |>
  unnest(pearson) |>
  mutate(cancer = fct_reorder(cancer, estimate),
         significance = case_when(p.value > .05 ~ 'NS',
                                  estimate > 0 ~ 'Positive',
                                  .default = 'Negative')) |>
  ggplot(aes(estimate, cancer)) +
  geom_col(aes(fill = significance)) +
  geom_text(aes(label = str_c('p=',signif(p.value, 2))), size = 2,
            x = .75, hjust = 'left') +
  expand_limits(x = 1.6) +
  geom_vline(xintercept = 0) +
  theme_jpub() +
  scale_fill_manual(values = c('grey','tomato')) +
  labs(title = 'TRPM2 correlation with inflammatory score in TME macrophages',
       x = 'Pearson correlation R value', subtitle = 'Zhang, et al. 2021. Cell')

publish_source_plot('pancancer.macro.m2.inflam.score.corr')

pan_cancer_m2inf_mf |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'KIDNEY', features.plot != 'inflam1') |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  labs(title = 'Kidney cancer macrophage')

publish_source_plot('kidney_cancer_macro_m2_dotplot')

pan_cancer_m2inf_mf |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'KIDNEY') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .2, 'TRPM2-hi MF', 'TRPM2-lo MF')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'TME macrophage clusters in kidney tumor',
       y = 'Imflammatory score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('kidney.cancer.mf.trpm2.inflam.correlation', width = 70)

### TRPM2 x cytokine ------------
cytk_genes <- c('TNF','IL6','IL1B','IFNB1','IL12A')

panc_mf_m2_cytk <- sobj_lst |>
  map(\(x)corr_in_celltype(x, 'Macro', c('TRPM2', cytk_genes)), .progress = T)

corr_genes_m2 <- function(df, x){
  df |>
    list_rbind(names_to = 'cancer') |>
    pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = c(cancer, id)) |>
    filter(!is.na(.data[[x]])) |>
    summarise(pearson = cor.test(TRPM2, .data[[x]]) |> broom::tidy(), .by = cancer) |>
    unnest(pearson) |>
    select(cancer, estimate, p.value)
  }

cytk_m2_macro <- cytk_genes |>
  set_names() |>
  map(\(x)corr_genes_m2(panc_mf_m2_cytk, x)) |>
  list_rbind(names_to = 'gene')

cytk_m2_macro |>
  ggplot(aes(gene, cancer, fill = estimate, size = -log10(p.value))) +
  geom_point(shape = 21) +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(cytk_m2_macro$estimate)) +
  labs(fill = 'Pearson correlation') +
  theme_bw()

cytk_m2_macro |>
  ggplot(aes(gene, cancer, fill = estimate)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(cytk_m2_macro$estimate)) +
  labs(fill = 'Pearson correlation',
       title = 'Correlation of TRPM2 with cytokine expression in TME macrophage') +
  theme_jpub(theme_classic)

publish_source_plot('pancancer.macro.cytokine.m2.corr.heatmap')

## monocytes ------------
pan_cancer_m2inf_mo <- sobj_lst |>
  map(\(x)corr_in_celltype(x, 'Mono', redo_pca = F), .progress = T)

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = c(cancer,id)) |>
  summarise(pearson = cor.test(TRPM2, inflam1) |> broom::tidy(), .by = cancer) |>
  unnest(pearson) |>
  mutate(cancer = fct_reorder(cancer, estimate),
         significance = case_when(p.value > .05 ~ 'NS',
                                   estimate > 0 ~ 'Positive',
                                   .default = 'Negative')) |>
  ggplot(aes(estimate, cancer)) +
  geom_col(aes(fill = significance)) +
  geom_text(aes(label = str_c('p=',signif(p.value, 2))), x = .1, color = 'white',
            size = 2, hjust = 'left') +
  geom_vline(xintercept = 0) +
  theme_jpub() +
  scale_fill_manual(values = c('grey','tomato')) +
  labs(title = 'TRPM2 correlation with inflammatory score in TME monocytes',
       x = 'Pearson correlation R value', subtitle = 'Zhang, et al. 2021. Cell')

publish_source_plot('pancancer.mono.m2.inflam.score.corr')

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'THCA', features.plot == 'TRPM2') |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  labs(title = 'THCA monocytes')

publish_source_plot('THCA_mono_m2_dotplot')

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'THCA') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .12, 'TRPM2-hi Mo', 'TRPM2-lo Mo')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'TME monocyte clusters in THCA',
       y = 'Imflammatory score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('THCA.mo.trpm2.inflam.correlation', width = 70)

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'KIDNEY', features.plot == 'TRPM2') |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  labs(title = 'Kidney cancer monocytes') 

publish_source_plot('kidney_cancer_mono_m2_dotplot')

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'KIDNEY') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .05, 'TRPM2-hi Mo', 'TRPM2-lo Mo')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(title = 'TME monocyte clusters in kidney cancer',
       y = 'Imflammatory score', x = 'Average expression of TRPM2') +
  theme_jpub

publish_source_plot('Kidney.cancer.mo.trpm2.inflam.correlation', width = 70)

#### OV-FTC ----------
pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'OV-FTC', features.plot == 'TRPM2') |>
  BubblePlot() +
  theme_jpub() +
  labs(title = 'Ovarian/Fallopian tube cancer monocytes') 

publish_source_plot('ovarian_cancer_mono_m2_dotplot')

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'OV-FTC') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .05, 'TRPM2-hi Mo', 'TRPM2-lo Mo')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'TME monocyte clusters in ovarian/fallopian tube cancer',
       y = 'Imflammatory score', x = 'Average expression of TRPM2') +
  theme_jpub()

publish_source_plot('ovarian.cancer.mo.trpm2.inflam.correlation', width = 70)

#### UCEC -------------
pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'UCEC', features.plot == 'TRPM2') |>
  BubblePlot() +
  theme_jpub() +
  labs(title = 'UCEC monocytes') 

publish_source_plot('uterine_cancer_mono_m2_dotplot')

pan_cancer_m2inf_mo |>
  list_rbind(names_to = 'cancer') |>
  filter(cancer == 'UCEC') |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = ifelse(TRPM2 > .05, 'TRPM2-hi Mo', 'TRPM2-lo Mo')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2, box.padding = .1) +
  stat_cor(size = 2) +
  labs(title = 'TME monocyte clusters in UCEC',
       y = 'Imflammatory score', x = 'Average expression of TRPM2') +
  theme_jpub()

publish_source_plot('uterine.cancer.mo.trpm2.inflam.correlation', width = 70)

### TRPM2 x cytokine ------------
cytk_genes <- c('TNF','IL6','IL1B')

panc_mo_m2_cytk <- sobj_lst |>
  map(\(x)corr_in_celltype(x, 'Mono', c('TRPM2', cytk_genes)), .progress = T)

cytk_m2_mono <- cytk_genes |>
  set_names() |>
  map(\(x)corr_genes_m2(panc_mo_m2_cytk, x)) |>
  list_rbind(names_to = 'gene')

cytk_m2_mono |>
  ggplot(aes(gene, cancer, fill = estimate)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(cytk_m2_mono$estimate)) +
  labs(fill = 'Pearson correlation',
       title = 'Correlation of TRPM2 with cytokine expression in TME monocyte') +
  theme_jpub(theme_classic)

publish_source_plot('pancancer.mono.cytokine.m2.corr.heatmap')

## kidney cancer -------------
sobj_lst |> names()

sobj <- sobj_lst[[2]]

sobj <- sobj |>
  mutate(manual_main = str_extract(MajorCluster, 'Mast|pDC|cDC|Mono|Macro'))

sobj |> DimPlot(group.by = 'manual_main', cols = 'Paired') +
  ggtitle('Kidney cancer') +
  theme_jpub(theme_classic)

publish_pdf('kidney.cancer.myeloid.umap.pdf', width = 60)

g1 <- sobj |>
  DotPlot(c('TRPM2','inflam1'), group.by = 'manual_main')

g1$data |>
  mutate(features.plot = ifelse(features.plot == 'inflam1',
                                'Inflammatory score', 'TRPM2')) |>
  BubblePlot() +
  theme_jpub() +
  labs(x = '', title = 'Kidney cancer TME') +
  RotatedAxis()

publish_source_plot('kidney.cancer.myeloid.m2.inflam.dotplot')

sobj |>
  DotPlot(c('TRPM2', cytk_genes), group.by = 'manual_main', cols = 'RdBu') |>
  pluck('data') |>
  BubblePlot() +
  theme_jpub() +
  labs(title = 'Kidney cancer TME') +
  RotatedAxis()

publish_source_plot('kidney.cancer.myeloid.m2.cytokine.dotplot')

sobj$tissue |> table()

sobj |>
  FeaturePlot('TRPM2', split.by = 'tissue', order = T,
              cols = c('lightgrey','tomato'))

sobj |>
  FindMarkersAcrossVar(split.by = 'manual_main', group.by = 'tissue',
                       ident.1 = 'T', features = 'TRPM2')

major_m2_bysample <- sobj |>
  DotPlot2d('TRPM2', group.x = library_id, group.y = MajorCluster) |>
  pluck('data')

major_m2_bysample |>
  mutate(group = str_extract(group.x, 'N$|T$')) |>
  ggplot(aes(group, avg.exp)) +
  geom_boxplot() +
  facet_wrap(~group.y, scales = 'free_y') +
  stat_compare_means(method = 't.test')

### M2h vs M2l ----------
#### macro -----------
sobj_mf <- sobj |>
  filter(str_detect(MajorCluster, 'Macro'))

sobj_mf <- sobj_mf |>
  quick_process_seurat(skip_norm = T, batch = 'library_id')

g1 <- last_plot()

g1 + theme_jpub(theme_classic) +
  ggtitle('Kidney cancer macrophages')

sobj_mf |> DotPlot(c('TRPM2','inflam1'), cols = 'RdBu')

sobj_mf <- sobj_mf |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(8,6), 'TRPM2-lo', 'TRPM2-hi'))

sobj_mf |>
  DotPlot(modc_marker, cols = 'RdBu')

sobj_mf |>
  bill.violin(inflam15$gene, trpm2_type, facet.ncol = 5, pubsize = T) +
  labs(x = '', fill = 'Subtype')

publish_pdf('kidney.cancer.mf.inflam.gene.violin.pdf', width = 90)

#### mono -----------
sobj_mo <- sobj |>
  filter(str_detect(MajorCluster, 'Mono'))

sobj_mo <- sobj_mo |>
  FindClusters(algorithm = 4, random.seed = 1)

g1 <- last_plot()

g1 + theme_jpub(theme_classic) +
  ggtitle('Kidney cancer macrophages')

sobj_mo |> DotPlot(c('TRPM2','inflam1'), cols = 'RdBu')

sobj_mo <- sobj_mo |>
  mutate(trpm2_type = ifelse(seurat_clusters %in% c(1,3,5,7), 'TRPM2-lo', 'TRPM2-hi'))

sobj_mo |>
  bill.violin(inflam15$gene, trpm2_type, facet.ncol = 5, pubsize = T) +
  labs(x = '', fill = 'Subtype')

publish_pdf('kidney.cancer.mo.inflam.gene.violin.pdf', width = 90)

mo_m2hvl_deg <- sobj_mo |>
  FindMarkers(group.by = 'trpm2_type', ident.1 = 'TRPM2-hi') |>
  as_tibble(rownames = 'gene')

mo_m2hvl_gsego <- mo_m2hvl_deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(T) |>
  clusterProfiler::gseGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', eps = 0)

mo_m2hvl_gsego <- mo_m2hvl_gsego |>
  clusterProfiler::simplify()

mo_m2hvl_gsego |>
  plot_enrichment() +
  theme_jpub() +
  labs(title = 'GO BP pathway enrichment of kidney cancer mono:\nTRPM2-high vs TRPM2-low')

publish_source_plot('kidney.cancer.mo.m2hvl.gsea', width = 80)

## MF+Mo ------------
pan_cancer_m2inflam <- pan_cancer_m2inflam |>
  set_names(pname) |>
  list_rbind(names_to = 'cancer')

pan_cancer_m2inflam |>
  write_source_csv('pan_cancer_trpm2_inflam_score')

pan_cancer_m2inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = c(cancer,id)) |>
  summarise(pearson = cor.test(TRPM2, inflam1) |> broom::tidy(), .by = cancer) |>
  unnest(pearson) |>
  mutate(cancer = fct_reorder(cancer, estimate)) |>
  ggplot(aes(estimate, cancer)) +
  geom_col() +
  geom_text(aes(label = str_c('p=',signif(p.value, 2))), x = .4) +
  expand_limits(x = .4) +
  geom_vline(xintercept = 0) +
  theme_bw() +
  labs(title = 'TRPM2 correlation with inflammatory score in TME macrophages/monocytes',
       x = 'Pearson correlation R value', subtitle = 'Zhang, et al. 2021. Cell')

pan_cancer_m2inflam |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = c(cancer,id)) |>
  filter(!str_detect(id, 'DC')) |>
  summarise(pearson = cor.test(TRPM2, inflam1) |> broom::tidy(), .by = cancer) |>
  unnest(pearson) |>
  ggplot(aes(estimate, cancer)) +
  geom_point() +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high)) +
  expand_limits(x = c(-1,1)) +
  geom_vline(xintercept = 0) +
  labs(title = 'TRPM2 correlation with inflammatory score in TME Macro/mono')



# convert logtpm back to counts --------------
esca_dense <- esca_dm |>
  transpose(keep.names = 'gene', make.names = 'index')

esca_meta |> filter(barcode == 'AAACGGGAGATCCTGT')

expm1(esca_dense$`AAACGGGAGATCCTGT-5`) / 7303 * 9676 |>
  tail(100)

median(esca_meta$n_counts)

esca_sparse <- esca_dense |>
  column_to_rownames('gene') |>
  as.sparse()
  
esca_7676 <- expm1(esca_sparse) / 7303

cell_count_sum <- esca_meta$n_counts

esca_nose <- esca_7676 |> nose()

esca_nose %*% Matrix::Diagonal(x = cell_count_sum) |> round()

esca_7676 |> lobstr::obj_size()

esca_7676 |> glimpse()

esca_round <- esca_7676 %*% Matrix::Diagonal(x = cell_count_sum)

esca_round |> nose()

esca_round@x <- esca_round@x |> round()

esca_round |> glimpse()

esca_round |> lobstr::obj_size()

colnames(esca_round) <- colnames(esca_7676)

esca_round |> DropletUtils::write10xCounts(path = 'mission/SLE_TRPM2_MfMo/data/z2021_pan_myl/esca.h5')

sobj <- esca_round |>
  CreateSeuratObject(meta.data = esca_meta)

h5_path <- list.files(pattern = 'h5')

h5_final <- pexpr |>
  str_replace('normalized.+', 'count.h5')

file.copy(h5_path, h5_final)
